rm(list=ls())

library(foreign)
library(locfit)
library(Hmisc)
setwd("/home/users/nycdoe.OU/functions")
source("rdrobust.R")
source("rdbwselect.R")
source("functions.R")
source("rdfunctions_ll_shu.R")
source("DRDplot.R")

setwd("/home/users/nycdoe.OU/results")
nyc <- read.dta(file="reg_v6.dta")
newdata<- nyc[c("std_mark","dist_cutoff","group","appyr","Regent_dbn","courseid")]
newdata$missing<-is.na(newdata$dist_cutoff)
newdata <-subset(newdata,missing==0)
smooth=4.5
table(newdata$appyr)
numd=2
#################################################################
for (d in 1:numd){
	mrdresults=matrix(NA,ncol=4,nrow=24)
    drdresults=matrix(NA,ncol=8,nrow=24)
	if (d==1){newdata2 <-subset(newdata,appyr>=2003&appyr<=2006)} else # EMA sample
	if (d==2){newdata2 <-subset(newdata,appyr>=2002&appyr<=2009)} # all years
	for (csid in 1:6){
		for (schid in 1:4){
		print(c(schid,csid))
		if (schid==1){
		data<-subset(newdata2,round(courseid)==csid)		
		} else {
		data<-subset(newdata2,round(courseid)==csid &round(group)==schid-1)	
		}	
ctff<-0
Xeval <- seq(-5, 5, length=51)  

################################################################
daty=data$std_mark
datx=data$dist_cutoff
NT=sum(datx>0)
NC=sum(datx<0)
ygrid=quantile(daty[abs(datx)<2],seq(0.001,0.999,length.out=999))
## Mean Change in Outcome
if (schid==1){
mrdbw.out = rdbwselect(daty,datx, bwselect = "IK")
bw_mrd=mrdbw.out$bws[1,1]*(NC+NT)^(1/5-1/smooth)
}
print(sum(abs(datx)<=bw_mrd))
results<-MRD_sharp(daty,datx,bw_mrd,Xeval)
mrdresults[(schid-1)*6+csid,]=c(bw_mrd,results$estimate,results$tstat,2-2*pnorm(abs(results$tstat)))

## DRD test
results<-DRDtest(daty,datx,ygrid,h=bw_mrd,bwsel="IK",undersmooth=smooth)
bw_drd=results$bw_iter
print(sum(abs(datx)<=bw_drd))
drdresults[(schid-1)*6+csid,]=c(results$bw,min(results$ksstat),max(results$ksstat),DRDpvalue(max(abs(results$ksstat))),results$bw_iter,min(results$ksstat_iter),max(results$ksstat_iter),DRDpvalue(max(abs(results$ksstat_iter))))

save.image(file=paste("nyc_results_sch",schid,"_csid_",csid,"_d_",d,".RData",sep=""))

## DRD plot
xlab="standardized Regents score"
ylab1="Conditional CDFs"
ylab2="Rescaled Difference in CCDFs"
if (csid==1){
main="Regents Math Score"	
} else if (csid==2){
main="Regents Advanced Math Score"	
} else if (csid==3){
main="Regents English Score"	
} else if (csid==4){
main="Regents Global History Score"	
} else if (csid==5){
main="Regents U.S. History Score"	
} else if (csid==6){
main="Regents Living Environment Score"	
} 
	
labels<-data.frame(xlab,ylab1,ylab2,main)
pdf(paste("/home/users/nycdoe.OU/results/graphs/nyc_results_sch",schid,"_csid_",csid,"_d_",d,".pdf",sep=""),width=7,height=5)
DRDplot(ygrid,results,labels,iter=1,type=2)
dev.off()
	}
}
print(mrdresults)
print(drdresults)
print(sort(mrdresults[7:24,4]))
print(sort(drdresults[7:24,8]))
}


